// This sampler uses trilinear interpolation
// The sample point should be in voxel space
template<typename T>
__device__ T sample(int nx, int ny, int nz, T x, T y, T z, T *f)
{
	T xf = min(max(x, T(0)), T(nx - 1));
	T yf = min(max(y, T(0)), T(ny - 1));
	T zf = min(max(z, T(0)), T(nz - 1));
	int xi = floor(xf);
	int yi = floor(yf);
	int zi = floor(zf);
	xf -= xi;
	yf -= yi;
	zf -= zi;
	if(xi == nx - 1) { --xi; xf += T(1); }
	if(yi == ny - 1) { --yi; yf += T(1); }
	if(zi == nz - 1) { --zi; zf += T(1); }
	int i0 = xi + yi * nx + zi * nx * ny;
	int i1 = i0 + 1;
	int i2 = i0 + nx;
	int i3 = i2 + 1;
	int i4 = i0 + nx * ny;
	int i5 = i4 + 1;
	int i6 = i4 + nx;
	int i7 = i6 + 1;
	T f01 = f[i0] * (T(1) - xf) + f[i1] * xf;
	T f23 = f[i2] * (T(1) - xf) + f[i3] * xf;
	T f45 = f[i4] * (T(1) - xf) + f[i5] * xf;
	T f67 = f[i6] * (T(1) - xf) + f[i7] * xf;
	T f0123 = f01 * (T(1) - yf) + f23 * yf;
	T f4567 = f45 * (T(1) - yf) + f67 * yf;
	return f0123 * (T(1) - zf) + f4567 * zf;
}
